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The goal of this investigation is to further develop nonlinear modal numerical simulation 
methods for prediction of geometrically nonlinear response due to combined thermal- 
acoustic loadings. As with any such method, the accuracy of the solution is dictated by the 
selection of the modal basis, through which the nonlinear modal stiffness is determined. In 
this study, a suite of available bases are considered including (i) bending modes only; (ii) 
coupled bending and companion modes; (iii) uncoupled bending and companion modes; and 
(iv) bending and membrane modes. Comparison of these solutions with numerical 
simulation in physical degrees-of-freedom indicates that inclusion of any membrane mode 
variants (ii - iv) in the basis affects the bending displacement and stress response 
predictions. The most significant effect is on the membrane displacement, where it is shown 
that only the type (iv) basis accurately predicts its behavior. Results are presented for beam 
and plate structures in the thermally pre-buckled regime. 
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Nomenclature 

mass matrix (physical and modal coordinates) 

proportional damping matrix (physical and modal coordinates) 
displacement response vector (physical and modal coordinates) 

force excitation vector (physical and modal coordinates) 
nonlinear restoring force (physical and modal coordinates) 

modal basis function matrix 
time 

linear, quadratic nonlinear, and cubic nonlinear modal stiffness coefficients 
identity matrix 
viscous damping factor 
undamped natural frequencies 


I. Introduction 

Future aircraft and launch vehicle structures will be exposed to increasing thermal and acoustic loads. In order 
to avoid the premature failure due to high cycle fatigue, fast and accurate analysis methods are required for use in 
the design process. The goal of this investigation is to further develop fmite-element-based nonlinear modal 
numerical simulation methods for prediction of geometrically nonlinear response due to combined thermal-acoustic 
loadings. These developments will ultimately form the basis of a computationally efficient thermal-acoustic fatigue 
design tool. 
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Reduced order analysis methods gain their computational advantage over direct numerical simulation by 
transforming the equations of motion in physical degrees-of-freedom to modal degrees-of- freedom. Hence system 
size may be significantly reduced. In the recent past, a significant amount of research has been performed in the 
area of reduced order methods development. These may be viewed as being in one of two categories; those in which 
the nonlinear modal stiffness is directly evaluated from the nonlinear finite element stiffness matrix (so-called direct 
methods), and those in which the nonlinear modal stiffness is indirectly evaluated. Direct methods are typically 
implemented in special purpose finite element codes in which the nonlinear stiffness is known. The work by Mei et 
al 1 for flat configurations and the expansion of that approach for curved panels by Przekop et al 2,3 are good examples 
of the direct method. Indirect stiffness evaluation methods are typically implemented for use with commercial finite 
element codes in which the nonlinear stiffness is unavailable. Examples of indirect stiffness evaluation approaches 
may be found in the work of McEwan et al 4 and Muravyov and Rizzi. 5 The latter implementation, as embodied in 
the code ELSTEP 6 , forms the basis for the current work. ELSTEP uses MSC.NASTRAN 2001 to evaluate the 
nonlinear stiffness. 

For both direct and indirect stiffness evaluation approaches, the crux of the problem lies in the selection of the 
proper basis, through which the nonlinear modal stiffness may be determined. If an insufficient basis is selected, 
then the predicted dynamic response of the reduced order model may significantly differ from that of the full 
nonlinear model. Thus, there is a need in all circumstances to assess the basis used through comparison of the 
predicted response with something other than a reduced order method, e.g. experimental data or numerical 
simulation in physical degrees-of-freedom. Recent work by Hollkamp et al 7 compared various reduced order 
approaches with experimental data. The use of experimental data, however, may be difficult, as both thermal and 
mechanical boundary conditions must be matched between the reduced order finite element analysis and experiment. 
Numerical simulation in physical degrees-of-freedom eliminates this difficulty as thermal and mechanical boundary 
conditions may be identically specified, though at the cost of a computationally intensive analysis. This approach 
has been successfully employed to validate a reduced order equivalent linearization analysis in the absence of 
thermal loading. 8 

The objective of the present work is to assess the effect of basis selection on thermal-acoustic random response 
prediction in the pre-buckled regime using a nonlinear modal simulation. For this class of problems, geometric 
nonlinearity couples the bending and membrane response, indicating a need to represent the effect of both in the 
basis selection. In an effort to determine which selection is most accurate, a suite of bases are considered including 
(i) bending modes only; (ii) coupled bending and companion modes; (iii) uncoupled bending and companion modes; 
and (iv) bending and membrane modes. Results from nonlinear modal simulation are compared with results 
obtained through numerical simulation in physical coordinates. Planar aluminum beam and plate structures are 
considered. 


II. Nonlinear Modal Simulation 

The nonlinear modal simulation analysis consists of several parts. The linear system matrices are first obtained 
from a commercial finite-element program. One or more methods, to be discussed, are next used to obtain a modal 
basis. Following a transformation of the nonlinear system to modal coordinates, the modal stiffness coefficients are 
evaluated and the resulting coupled system of equations is numerically integrated to obtain the modal displacement 
time history. These are transformed back to physical coordinates for post-processing, including stress recovery. 

A. Modal Coordinate Transformation 

The equations of motion of the nonlinear system in physical degrees-of-freedom may be written as 

MX(t) + CX(0 + F nl (A(0) = F(t) . (1) 

A set of coupled modal equations with reduced degrees-of-freedom is first obtained by applying the modal 
coordinate transformation X = to Eq. (1). The modal basis function matrix O is typically formed from the 

eigenvectors obtained from Eq. (1) using only the linear stiffness. For flat isotropic structures, these may include 
any combination of bending and membrane modes. In lieu of membrane modes, the modal basis may include 
“companions” related to the membrane response, as discussed in the next section. Generally, a small set of L basis 
functions are included giving 

Mq + Cq + F NL (q l ,q 2 ,...,q L ) = F. (2) 
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The time dependence of q , F and F NL is implied, and 

M = O r MO = |7 J C = O r CO = T 2 C r a> r J F NL = ® T F m # = O 1 /’ . (3) 

B. Modal Basis Calculation 

For the problems of interest in this paper, both bending and membrane behavior should be included in the basis 
selection since the large deflection nonlinearity couples their response. The basis functions may be determined via 
several methods. Basis functions corresponding to the bending and membrane response may be determined through 
solution of the linear eigenvalue problem. Other basis functions corresponding to the membrane response induced 
by bending-membrane coupling may be determined via alternative approaches. Lastly, while not the subject of this 
paper, basis functions may also be determined via experiment or hybrid scheme. 

1. Linear Eigenvectors 

In this study, only isotropic materials or laminated composite materials having zero laminate coupling stiffness 
([B]=0) are considered. Further, the current work is restricted to planar structures, i.e. flat beams and plates. The 
linear eigenvectors obtained from Eq. (1) (using only the linear stiffness) are therefore uncoupled and are either 
associated with low-frequency bending modes or high-frequency membrane modes. The selection of which bending 
modes to include depends on both the excitation bandwidth and the loading distribution. The identification of each 
selected bending mode is straightforward as they are typically the lowest frequency modes. The selection of which 
membrane modes to include is less apparent. All or most membrane modes will likely fall outside the excitation 
bandwidth. Furthermore, the identification of any particular membrane mode may be difficult and time consuming. 
Nevertheless, a reasonable starting point is to select the lowest membrane modes that are consistent with the loading 
distribution. Inclusion of both these eigenvectors in the modal basis, either independently or in pairs, is 
subsequently referred to as the bending and membrane mode (BM) case. Inclusion of only the bending eigenvectors 
will be referred to as the bending mode only (B) case. 

In this study, the mass-normalized linear eigenvectors without the effect of thermal stress were obtained using 
MSC.NASTRAN normal modes analysis (solution 103). These are the so-called “cold” modes. The mass- 
normalized linear eigenvectors with the effect of thermal stress, i.e. “hot” modes, requires MSC.NASTRAN 
nonlinear static analysis (solution 106). 

2. Companion Basis Functions 

An alternative approach to using membrane modes is the use of so-called companion 7 (or dual 9 ) modes. These 
modes represent the membrane behavior resulting from bending due to bending-membrane coupling. Previous 
authors utilized quasi-static approaches to determine the companion mode via nonlinear static analyses. For each 
companion mode, Hollkamp et al 7 prescribed a displacement field corresponding to the bending mode of interest. 
Mignolet et al 9 applied two loadings, each having the distribution of a particular bending mode shape but with 
different magnitudes, to obtain the companion for that bending mode. In either case, the resulting modes contained 
only the membrane related behavior, but not the bending. 

A new method is now presented for computing the companion mode using a dynamic analysis. An initial stress- 
free imperfection is first introduced in the flat structure to couple the bending-membrane response. The magnitude 
of the imperfection is chosen to be very small such that a normal modes analysis yields virtually the same bending 
eigenvalues and eigenvectors as that of the flat structure. An engineering rationale for specifying the shape of the 
imperfection is not evident. In this paper, the shape was chosen to be that of the first bending mode. It was found, 
however, that other imperfection shapes resulted in different companion parts of the eigenvectors. The bending and 
high frequency membrane components remained identical to that of the flat structure regardless of the shape of the 
imperfection. For the cold modes, MSC.NASTRAN solution 103 was again used to compute the mass-normalized 
eigenvectors, which now contain both the bending and membrane behaviors, but at the natural frequencies of the 
original flat structure. Direct inclusion of these eigenvectors in the modal basis is subsequently referred to as the 
coupled bending and companion mode (CBC) case. A more consistent usage, however, is to separate the degrees- 
of-freedom (DoFs) associated with the bending and membrane behaviors. In practice, since the bending behavior is 
unchanged, the original low-frequency bending modes are retained. The bending DoFs are set to zero in the newly 
obtained eigenvector to obtain the dynamic companion. Inclusion of both these eigenvectors in the modal basis, 
independently or in pairs, is subsequently referred to as the uncoupled bending and companion mode (UBC) case. 
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3. Comparison of Membrane and Companion Modes 

Consider a clamped-clamped aluminum beam measuring 18-in. x 1-in. x 0.09-in (/ x w x h) with the following 
material properties: 

E = 10.6x10 6 /7«, G = 4.0 xlO 6 psi, p = 2.588 x KT* . (4) 

The beam was modeled using seventy-two CBEAM elements. The first four symmetric bending modes 
(eigenvectors 1, 3, 7 and 10) are at natural frequencies of 58, 312, 770 and 1430 Hz, respectively. The first four 
anti-symmetric membrane modes (eigenvectors 47, 86, 127, 140) are at natural frequencies of 11.2, 22.5, 33.6 and 
44.8 kHz, respectively. These are plotted in Figure 1 - Figure 4. Also shown are the static companion modes 
obtained by Mignolef s approach 9 and the dynamic companion modes obtained by the method outlined above. It is 
clear that both static and dynamic companion modes significantly differ in shape relative to the membrane mode. 
The effect of their inclusion in the modal basis on the stiffness coefficients is next considered. 

C. Indirect Stiffness Evaluation Method 

The indirect stiffness evaluation method previously developed without the effect of thermal loading 5 was further 
developed for thermal loadings. 9 To summarize, the nonlinear force vector in Eq. (2) may be written in the form 

L L L L L, Z, 

F NL (q t ,q 2 ,...,q L ) = Y J d[q. + Z Z A + Z Z Z r = l, 2,..., L (5) 

j = 1 j= 1 k=j j = 1 k=j l=k 

reducing the problem of determining the nonlinear stiffness from one in which a large set of simultaneous nonlinear 
equations must be solved to one involving simple algebraic relations. For the prescribed static displacement fields 9 

V, = +< H V-2 = -<tw, x ci = +< l>,?, ( 6 ) 

the nonlinear forces evaluated using the MSC.NASTRAN nonlinear static solution (solution 106) are given as 

E, = (+H ) = + [< ] + [<> ] C k C k + Em ] 

(-H ) = - [< ] ?1 + [<> ] ?1?1 - [*m ] ?.?.?. • O) 

E, = (+4tf ,) = + [< ] i + [a n ] ii - \b[ u ] m 

The first two of Eqs. (7) are used to evaluate | ~a r n J . The first and the third of Eqs. (7) are then used to determine 

\jd[ J and \jbf ] . It should be noted that the specified modal displacements in Eq. (6) are single scalar quantities 
since the analysis is static. These modal displacements are specified such that the magnitude of the prescribed 
physical displacement field X is physically meaningful. It has been shown that the stiffness coefficients obtained 

via this approach are not sensitive to the particular value of modal displacement specified. 5 The last of Eqs. (6) are 
obtained by specifying a scaled modal displacement such that q - aq . In this paper, a value of a =1.25 is used. 

This procedure for evaluating the linear stiffness coefficients [^ J is advantageous over the earlier 

implementation 5 using the MSC.NASTRAN linear static solution (solution 101), as the linear solution does not 
include in-plane displacements, and therefore is not capable of handling thermal loads. The remaining nonlinear 
coefficients are evaluated through an additional series of nonlinear static solutions. 5 

1. Comparison of Modal Stiffness Coefficients 

For the clamped-clamped beam structure previously considered, the effect of including a particular set of modes 
in the modal basis is shown in Table 1-3. Table 1 shows the diagonal terms of the linear stiffness coefficients for a 

set consisting of the first symmetric bending modes ( d\ - d\ ) plus up to four membrane or companion modes ( d\ - 
d\ ). In each case, the first four linear stiffness coefficients are equal to the eigenvalues of the selected bending 
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modes. For the BM case, the second four linear stiffness coefficients correspond to the eigenvalues of the selected 
membrane modes. The CBC case lacks the second set of four coefficients as only four modes are included in this 
basis, i.e. four coupled bending and companion modes. Both versions of the UBC case indicate low frequency 
behavior for the second set of four (companion) modes. In all cases but CBC, the off-diagonal terms are effectively 
zero (not shown). 

Since the beam structure is flat, all quadratic stiffness coefficients are nominally zero and are not shown. Table 2 
shows the diagonal terms of the cubic stiffness coefficients for the same modal bases as above. In each case, all 
coefficients corresponding to the bending modes (1-4) are identical, indicating that the effect of the different modal 
bases is reflected in the off-diagonal terms that couple the various modes. This is more apparent in Table 3, which 
shows selected off-diagonal stiffness coefficients corresponding to the first (1 st symmetric bending) and fifth (1 st 
membrane or companion) modes. There is a much greater coupling of bending and membrane modes in the BM 

basis than any of the other bases, as shown by stiffness terms b' {55 and b s U5 . Additionally, the relationship between 
non-zero stiffness terms, 5 

bL=b 5 n5 ( 8 ) 

is approximately maintained for the BM and UBC bases, and somewhat less so for these particular terms in the 
static UBC basis. Further discussion regarding the choice of basis on the response prediction is reserved until Sec. 
III. 

D. Numerical Integration and Element Stress Recovery 

Having the nonlinear force vector in Eq. (5) fully defined, the coupled modal nonlinear equations of motion in 
Eq. (2) are numerically integrated using a fourth-order Runge-Kutta method. 10 The resulting modal displacement 
time histories are transformed back to physical coordinates using the inverse modal transformation. A fixed time 
integration step of between 50 jUs (for the lowest excitation levels) and 10 jus (for the highest excitation levels) was 
used. 

In contrast to mapping techniques relating stress to physical DoFs, 7 ’ 9 the approach taken here recovers the 
element stresses by post-processing the nodal physical DoFs directly within the finite element program. Since the 
stress post-processor of the finite element program is used, the stresses are identical to those that would have been 
obtained via a standard finite element analysis in physical DoFs. The appeal of direct stress recovery over mapping 
techniques is that it eliminates any issues concerning the accuracy of the mapping, and eliminates the need for 
estimating nonlinear mapping functions for each stress quantity of interest. The disadvantage is that it can be 
computationally costly. 

In principal, the approach taken is simple. For a single element stress recovery, a finite element mesh is made 
consisting of two elements. The stress recovery element is of the same type and has the same properties of the 
element of interest. For a particular output time step, the element physical DoFs, obtained via the nonlinear modal 
simulation method, are applied to each element node as prescribed displacement fields in the MSC.NASTRAN 
nonlinear static solution. A nuance of MSC.NASTRAN requires the attachment of a dummy element and load 6 to 
the stress recovery element, in order to calculate the resulting element nodal forces and recover the element stress. 
By repeating this operation for each output time step, the stress time history is determined. If stress recovery is 
required for a region or the entire structure, that portion of the original model is used and a dummy element is 
attached. The prescribed displacement field spans the entire region of interest. 

E. Convergence 

There are several aspects relating to eigenvalue and modal convergence that warrant discussion. A means of 
assessing the adequacy of the finite element mesh within a given frequency range is to check convergence of the 
eigenvalues. In the analyses performed herein, eigenvalue convergence was checked through the applied loading 
bandwidth. For the beam, natural frequencies of the membrane modes in the BM basis (11.2 - 44.8 kHz) were 
significantly higher than the loading bandwidth of 0-1500 Hz. Further, it is known that the membrane response has 
double the frequency content of the loading due to the fact that the membrane vibration undergoes two cycles for 
every one cycle of bending. In this study, no attempt was made to further refine the finite element mesh to achieve 
eigenvalue convergence to either double the maximum loading frequency or to the membrane mode eigenvalue. 
The implication of the latter is that membrane eigenvectors included in the basis will be less accurate than the 
bending eigenvectors. 

For any given modal basis variant, modal convergence (the effect of adding additional modes to the basis) must 
be assessed through convergence of the predicted response in either modal or physical coordinates. Since the 
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stiffness coefficients are independent of each other, the inclusion of additional modes in the basis only augments the 
size of the modal stiffness matrices. The addition of the n th basis does not affect the magnitudes of the coefficients 
evaluated using n-1 basis functions. 

III. Displacement and Stress Response Results 

A. Numerical Simulation in Physical Coordinates 

Numerical simulation analysis in physical coordinates served as the basis for comparison of results from 
nonlinear modal simulation analyses. The finite element program ABAQUS (version 6.4) was to generate nonlinear 
displacement and stress time histories. The double precision explicit integration scheme with an adaptive time 
integration step (referred to as “element by element” in ABAQUS) was utilized for all analyses. The mass 
proportional damping factor used was the same as that specified in the nonlinear modal simulation analyses. The 
choice of ABAQUS over alternative analysis methods, e.g. MSC.NASTRAN nonlinear transient solution (solution 
129), was made based on its superior ability to simulate very long response histories. This choice introduced an 
inconsistency of element formulations used between the nonlinear modal and physical simulations. However, nearly 
identical nonlinear static analysis results from MSC.NASTRAN and ABAQUS (not shown) established the 
consistency of displacement and stress. 

B. Load Generation, Stress Post-Processing, and Ensemble Averaging 

The same loading time history, ensemble averaging, and stress post-processing was used for the numerical 
simulation analysis in physical and modal coordinates. The loading time histories were generated by summing equal 
amplitude sine waves, each with random phase, within the specified bandwidth using a discrete inverse Fourier 
transform. This procedure was identical to that used in previous work 8 by the authors, so further details are omitted 
for brevity. The loading produced by this method has a Gaussian distribution. A sharp roll-off of the input 
spectrum practically eliminates excitation of the structure outside the frequency range of interest. 

The longitudinal stresses at the upper and lower surfaces were computed using the method previously described 
for the nonlinear modal simulation, or directly by ABAQUS for the simulation in physical coordinates. The bending 
and membrane components were separated by averaging the upper and lower surface stresses to obtain the 
membrane component, and then subtracting that from the surface stress to obtain the bending component. This 
separation allowed greater insight into the effect of the various modal bases on the stress response. 

A range of load levels was considered to span the response regime from nearly linear to highly nonlinear. For 
each load level, ten ensembles of displacement and stress response were generated lasting 2.1384s each. The first 
0.5s of each response record was discarded to remove the initial transient response, 8 resulting in response histories 
of 1.6384s in duration. In the case of the nonlinear modal simulation, the transient response was removed from the 
displacement time history prior to stress recovery. For each simulation, the displacement and stress response was 
stored at every 50 jus , giving time records of 32,768 points. A 32,768-point FFT was subsequently used to compute 
the power spectral density (PSD) function. 

C. Beam Results 

The clamped-clamped beam structure previously discussed was first considered. The MSC.NASTRAN model 
consisted of seventy-two CBEAM elements. Nonlinear modal simulations were run using four of the five modal 
bases from Table 1 - 3; B, BM, CBC and UBC. Simulations for the fifth basis, UBC static, were run for a subset of 
conditions. The ABAQUS model consisted of seventy-two B21 elements. Damping was chosen to be sufficiently 
high so that a good comparison could be made at the peaks of the PSD. A level of mass proportional damping was 
specified corresponding to critical damping of 2.0% for the first symmetric bending mode. The beam was subjected 
to a uniformly distributed transverse loading with a bandwidth of 1500 Hz. Two thermal conditions were 
investigated; a room temperature condition and an elevated temperature ( AT = 5°F ) below the critical thermal 

buckling temperature ( A T CR = 6.64 °F ). 

Displacement results are presented at the 23% span (4.125-in. from the clamped end), and element stresses at 
21% span (3.75-in. from the clamped end). These locations help to elucidate the benefits and liabilities of the 
various modal bases under severe conditions as the membrane stress component is more significant there relative to 
the bending stress component than at other locations along the length, e.g. the clamped end or mid- span. 
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1. Room Temperature 

The root-mean-square (RMS) transverse displacement is shown in Figure 5 for a range of loadings. Nonlinearity 
is seen to increase with increasing load. The best comparison with numerical simulation in physical DoFs is 
achieved using the BM basis. Both UBC and CBC bases are seen to offer some improvement relative to the bending 
only basis. The RMS membrane displacement is shown in Figure 6. Compared to the ABAQUS results, the level of 
agreement achieved using the BM basis is excellent. The largest contribution to the RMS membrane displacement 
using the UBC comes from the low frequencies of the companion modes (see Table 1). This will be further 
substantiated by the membrane displacement PSDs to follow. In an RMS sense, the membrane displacements 
obtained using the CBC basis were negligible. Note that the bending only basis is incapable of generating any 
membrane displacement. 

The transverse displacement PSDs are shown in Figure 7 for the lowest loading of 7.2x1 0" 3 lb/in. For this level, 
all modal bases give essentially the same results as the numerical simulation in physical DoFs. The membrane 
displacement PSDs for this load level, shown in Figure 8, offer a different perspective. The BM basis is the only 
one to compare well with the ABAQUS results. The CBC basis is inaccurate on two accounts; the magnitude is 
clearly incorrect and the shape of the frequency response mimics that of the bending response in Figure 7. It comes 
as no surprise that both UBC bases compare poorly. The magnitude of bending-companion coupling in the cubic 
stiffness coefficients are orders of magnitude lower than the bending-membrane coupling in the BM basis, see Table 
3, and therefore account for the near negligible displacement response other than at very low frequency. Since the 
two UBC bases give very similar results, differences between the stiffness terms in Table 3 are largely irrelevant and 
only affect some relative peak values and not the overall behavior. Interestingly the character of the peaks is 
captured in a qualitative sense. Though the loading spectrum falls off sharply at 1500 Hz, the nonlinear modal 
simulations of membrane displacement response using the B, BM, and UBC bases extend to twice that frequency 
because of period doubling. Membrane displacement response computed using the CBC basis does not exhibit 
period doubling. 

The bending stress PSDs at the lowest excitation level are shown in Figure 9. All modal bases compare 
favorably with the results from numerical simulation in physical DoFs at this level. The membrane stress PSDs are 
shown in Figure 10. With regard to the bending only modal basis, the membrane stress arises from the bending- 
membrane coupling in the stress recovery calculation, even though a zero membrane displacement is imposed. 
Results computed using the bending only and UBC bases over-predict the membrane stress at the first period- 
doubled peak, and capture the behavior up to the third peak at about 370 Hz. Above that frequency, the magnitudes 
differ significantly from the ABAQUS results. The same is true of the CBC basis, except that it additionally 
exhibits peaks at the bending frequencies. The B, UBC and CBC membrane displacements are vastly different than 
that of the BM basis, yet the membrane stresses have similar characteristics. It thus appears that the membrane 
stress in the low frequency range is more dependent upon the bending displacement than on the membrane 
displacement. The BM results compare favorably to a frequency of about 700 Hz, above which more significant 
differences are noted. This suggests that additional membrane modes in the BM basis are required to achieve 
membrane stress comparisons on a par with the membrane displacement comparisons. Above approximately 1150 
Hz, results from all modal bases give similar results, and differ significantly with those from numerical simulation in 
physical DoFs. 

The following results and subsequent figures for the high loadings have been smoothed to more clearly observe 
the behavior. At the highest loading of 0.9216 lb/in., the transverse displacement PSD exhibits the peak spreading 
and shifting characteristic of spring hardening nonlinearity, see Figure 11. For this case, the bending only basis 
indicates a higher degree of nonlinearity than the ABAQUS results, as evidenced by the first and second peaks 
shifted to higher frequencies, greater broadening, and smaller magnitudes. The UBC and CBC results are similar to 
the bending modes only, except that the fundamental mode is more accurately captured. The BM basis results 
compare well across the frequency range, with somewhat broader peaks than indicated by the numerical simulation 
in physical DoFs. The membrane behavior shown in Figure 12 follows the observations made at the lower level. 
The excellent agreement between BM basis and ABAQUS results demonstrate the efficacy of the BM basis in the 
highly nonlinear regime. 

Bending and membrane stress PSDs are shown for the highest loading in Figure 13 and Figure 14, respectively. 
The BM basis most closely compares with the ABAQUS bending stress across the frequency range. The bending 
only, UBC and CBC bases significantly under-predict the bending stress of the first and second modes, and above 
about 1100 Hz indicate different character than the BM basis and ABAQUS results. Further, the UBC and CBC 
bases exhibit a non-physical twin peak in the vicinity of the first mode. In Figure 14, the BM basis results compare 
favorably with the ABAQUS results across the frequency range. The bending only, UBC and to a lesser extent CBC 
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bases over-predict the membrane stress magnitude at the peak as well as the width of the peak. Above about 1100 
Hz, the character substantially differs from the ABAQUS results. 

2. Elevated Temperature 

The bending and membrane displacement and stress behavior for a beam under combined thermal and acoustic 
loading is next considered. The dimensions, material properties and boundary conditions used were the same as 
those previously used. A coefficient of thermal expansion of 1.24x1 0" 5 1/°F was additionally specified. The thermal 
loading imposed was such that the beam remained in the pre-buckled regime. 

The modes included in each modal basis were identical to those used for the room temperature case. Results 
computed using these “cold” modes compared more favorably with numerical simulation in physical DoFs than did 
those computed with “hot” modes. In the following, only results computed using the “cold” modes are presented. 

The RMS transverse displacements are shown in Figure 15. Unlike the room temperature case, transverse 
displacements for the elevated temperature case indicate no clear winner across the loading range. Depending on 
the level of dynamic loading, results from a particular basis may compare more or less favorably than the others. 
The membrane displacements shown in Figure 16, however, indicate the best comparison with the BM basis. Like 
the room temperature case, the largest contribution to the RMS membrane displacement using the UBC comes from 
very small eigenvalues of the companion modes. At the lowest dynamic loading, the trends observed in the 
transverse displacement PSD (Figure 17) follow those from the room temperature case. Similarly, the membrane 
displacement PSD (Figure 18) shows good agreement only between the BM basis and ABAQUS, with the most 
significant discrepancy in magnitude at the first peak. The bending stress PSD (Figure 19) and membrane stress 
PSD (Figure 20) have similar trends to those observed for the room temperature case. 

At the higher load level of 0.9216 lb/in., the transverse displacement PSDs from all modal bases capture the 
fundamental mode fairly well, as shown in Figure 21. Only the BM basis describes the second mode behavior well, 
with the other bases indicating reduced amplitude and greater broadening. The membrane displacement PSD 
computed using the BM basis is the only one to compare well with the ABAQUS results, see Figure 22. The 
bending and membrane stress PSDs are shown in Figure 23 and Figure 24, respectively. Like the room temperature 
case, the bending stress results obtained using the BM basis are the most accurate, with the bending modes, UBC 
and CBC bases under-predicting the peaks. Further, the UBC and CBC bases again exhibit the non-physical double 
peak at the first mode. The low frequency membrane behavior is most accurately captured using the bending only 
basis, followed by the BM basis, which indicates greater peak broadening. The BM basis most accurately predicts 
the high frequency membrane behavior. 

This set of results demonstrates the effect of modal basis selection on the nonlinear modal simulation of thermal- 
acoustic response in the pre-buckled regime. The quality of displacement and stress PSDs are consistent with the 
room temperature case. 

D. Plate Results 

A simply supported aluminum plate was next considered to help understand the complexities of modal basis 
selection for a two-dimensional structure. The plate measured 10-in. x 14-in. x 0.04-in (/ x w x h) with the material 
properties indicated in Eqs. (4) and was subjected to a uniformly distributed acoustic load with a bandwidth of 1024 
Hz. 

The MSC.NASTRAN model consisted of 560 CQUAD4 elements measuring 0.5 -in. x 0.5-in. The ABAQUS 
model was similarly constructed using S4R elements. Nonlinear modal simulations were run using three modal 
bases; B, BM and UBC. For the bending modes only basis, the first five doubly symmetric bending modes 
(eigenvectors 1, 4, 8, 11, 12) were included. These occurred at natural frequencies of 58, 217, 373, 529 and 543 Hz. 
For the BM basis, the first four anti- symmetric membrane modes (eigenvectors 276, 281, 364, 373) were 
additionally included for a total of 9 modes. These occurred at natural frequencies of 15.2, 15.5, 21.3 and 22.1 kHz. 
The selection of these particular membrane modes was somewhat arbitrary as there was no guiding rationale. Their 
identification amongst the hundreds of eigenvectors was also difficult. For the UBC basis, the five companions of 
the bending modes were included for a total of 10 modes. 

For this structure, an alternative nonlinear modal simulation using a direct stiffness evaluation method 2 was also 
available. This approach neglects in-plane inertia and statically condenses the membrane displacement into the 
bending displacement. Since this approach requires manipulating the nonlinear stiffness formulation, it is not 
possible to implement in general purpose commercial codes. To avoid confusion, this nonlinear modal simulation 
will subsequently be referred to as “static physical condensation” method. A quarter-plate finite element model was 
constructed from 140 0.5-in. square modules, each having two MIN3-derived triangular elements with first order 
shear deformation. 2 Symmetry conditions were imposed along the two centerlines. In the stress calculations to 
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follow, the element stresses of the two triangular elements within the module were averaged to obtain the stress at 
the center of the square module. This allowed the stress from the static physical condensation method to be 
compared with the nonlinear modal simulations and simulations in physical DoFs. 

Mass proportional damping corresponding to critical damping of 2.0% for the first bending mode was used in all 
three analyses. Displacement results are presented at the (%, %) point at the coordinate position (2.5-in., 3.5-in.), 
and element stresses nearest the (%, l A) point at the coordinate position (2.25-in., 3.25-in.). The stresses presented 
are along the short (10-in.) direction. 

The transverse displacement PSDs for an acoustic loading of 106 dB are shown in Figure 25. All nonlinear 
modal simulations including the static physical condensation method agree well with the ABAQUS results. 
Differences in the peak amplitudes and frequencies, and the lack of higher frequency modes in the modal simulation 
are clearly attributable to lack of converged meshes and to modal truncation. The corresponding membrane 
displacements PSDs are shown in Figure 26. Results computed using the BM basis compare very well with the 
numerical simulation in physical DoFs through about half the frequency range, above which additional membrane 
modes in the basis appear to be required. Results from the static physical condensation capture the first period- 
doubled peak, and attempt to capture the second and third period-doubled peaks. However, membrane 
displacements calculated using this method lack the essential character elsewhere. The UBC results indicate near 
zero membrane response, except at very low frequency. Figure 27 shows the bending stresses computed using each 
of the modal simulations compares favorably at the peaks, with differences elsewhere again attributable to lack of 
mesh convergence. Interestingly, the membrane stresses from the modal simulations agree fairly well with each 
other, and begin to differ substantially from the ABAQUS results above 300 Hz, see Figure 28. Since both the 
bending only and UBC modal bases produce membrane stresses on a par with those from the BM basis and static 
physical condensation, it further substantiates the contention that the membrane stress is more highly dependent on 
the bending displacement than on the membrane displacement. 

For an acoustic loading of 160 dB, a high degree of nonlinearity is observed in the transverse displacement PSDs 
shown in Figure 29. What was the first resonant peak at 58 Hz has substantially broadened and shifted to about 300 
Hz. Results from the BM basis and static physical condensation methods compare the best with numerical 
simulation in physical DoFs at the peak. The merits of basis selection or method used for the membrane 
displacement PSD is ambiguous, see Figure 30. At best, it can be seen that the UBC basis gives an effectively zero 
membrane response. The precise cause for this lack of agreement is presently unknown, although truncation of 
membrane modes in the BM basis is likely responsible for some of its difference. Finally, all of the modal 
simulations capture the peak bending stress response at about 300 Hz, as shown in Figure 31. The off-peak 
behavior, particularly at the higher frequencies, substantially differs from the numerical simulation in physical DoFs. 
The membrane stress PSDs from all modal simulations agree fairly well across the frequency range and capture the 
period-doubled peak at around 600 Hz, see Figure 32. 

In general, results obtained from any of the modal simulations of the plate response are of lesser quality than the 
modal simulations of the beam response, relative to their respective numerical simulations in physical DoFs. Some 
differences are clearly identifiable through observation of the low response level, including modal truncation and 
inadequate mesh convergence. These are likely to have also affected the quality of results at the higher levels, but in 
a less transparent way. In addition, since the membrane modes for the plate are more complicated than those of the 
beam, e.g. two-dimensional versus one-dimensional, the likelihood of selecting the proper modes for inclusion in the 
BM basis is lessened. Thus, the advantages of using the BM basis are diminished by practical considerations that 
may be difficult to overcome. 


IV. Conclusions 

The effect of modal basis selection on the displacement and stress response predicted using a nonlinear modal 
simulation was investigated. For the beam, it was found that 

• the bending and membrane modal basis offered the only accurate prediction of membrane displacement and 
the best prediction of bending and membrane stress, 

• the bending modes only basis offered the second best approach as it did not produce anomalous behavior in 
the displacement or stress PSDs, did not increase the size of the coupled set of modal equations, and did not 
have the added complication of identifying membrane modes or computing companion modes, 

• either variant of the companion modes basis (CBC or UBC) produced non-physical behavior in the 
displacement and stress PSDs, and therefore their use is not recommended, 
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• there is no advantage for computing companion modes dynamically (as presented in this paper) versus 
statically , 7,9 nor is there an engineering rationale for selection of eigenvector (for dynamic) or prescribed 
displacement or loading (for static), 

• there is no established criteria for selection of which modes to include in the basis for any of the membrane 
variants (BM, UBC, CBC). 

Results for the beam demonstrated that it is possible to construct an effective bending and membrane modal basis 
through the use of engineering judgment and some effort to identify the particular membrane modes selected. For 
the plate, it was additionally found that 

• the selection and identification of modes was more difficult and the accuracy of the predictions suffered, 

• the static physical condensation approach was able to accurately predict the first mode of the membrane 
displacement. It has the advantages of not increasing the system size and not having to identify or calculate 
membrane or companion modes, but is limited to implementation in special purpose codes. 

In general, increasing structural complexity requires additional bending and membrane modes in the basis to obtain 
the same level of accuracy. 

Finally, for the problems considered (those with nonlinear bending-membrane coupling), the calculation of 
membrane stress was more dependent upon the bending displacement than on the membrane displacement, 
particularly at low frequencies. Two areas for future study include the effect of basis selection on fatigue life 
predictions, and the effect of basis selection for problems that exhibit other types of nonlinearity, e.g. large rotation 
for cantilevered structures. 


References 

! Mei, C., Dhainaut, J.M., Duan, B., Spottswood, S.M., and Wolfe, H.F., "Nonlinear random response of composite panels 
in an elevated thermal environment," Air Force Research Laboratory AFRL-VA-WP-TR-2000-3049, Wright-Patterson Air 
Force Base, OH, October 2000. 

2 Przekop, A., "Nonlinear response and fatigue estimation of aerospace curved surface panels to combined acoustic and 
thermal loads," Ph.D. Dissertation, Old Dominion University, Norfolk, VA, 2003. 

3 Przekop, A., Guo, X., Azzouz, M.S., and Mei, C., "Reinvestigation of nonlinear random response of shallow shells using 
finite element modal formulation," Proceedings of the 45th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics 
and Materials Conference , AIAA-2004-1553, Palm Springs, CA, 2004. 

4 McEwan, M.I., Wright, J.R., Cooper, J.E., and Leung, Y.T., "A finite element/modal technique for nonlinear plate and 
stiffened panel response prediction," Proceedings of the 42nd AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, 
and Materials Conference , AIAA-2001-1595, Seattle, WA, 2001. 

5 Muravyov, A. A. and Rizzi, S.A., "Determination of nonlinear stiffness with application to random vibration of 
geometrically nonlinear structures," Computers and Structures , Vol. 81, No. 15, 2003, pp. 1513-1523. 

6 Rizzi, S.A. and Muravyov, A. A., "Improved equivalent linearization implementations using nonlinear stiffness 
evaluation," NASA TM-2001-210838, March 2001. 

7 Hollkamp, J.J., Gordon, R.W., and Spottswood, S.M., "Nonlinear sonic fatigue response prediction from finite element 
modal models: a comparison with experiments," Proceedings of the 44th AIAA/ASME/ASCE/AHS/ASC Structures, Structural 
Dynamics, and Materials Conference , AIAA-2003-1709, Norfolk, VA, 2003. 

8 Rizzi, S.A. and Muravyov, A.A., "Comparison of nonlinear random response using equivalent linearization and numerical 
simulation," Structural Dynamics: Recent Advances, Proceedings of the 7th International Conference , Vol. 2, Southampton, 
UK, 2000, N.S. Ferguson, et al. (ed.), pp. 833-846. 

9 Mignolet, M.P., Radu, A.G., and Gao, X., "Validation of reduced order modeling for the prediction of the response and 
fatigue life of panels subjected to thermo-acoustic effects," Structural Dynamics: Recent Advances, Proceedings of the 8th 
International Conference , Southampton, UK, 2003, M.J. Brennan, et al. (ed.). 

10 Press, W.H., Teukolsky, S.A., Vetterling, W.T., and Flannery, B.P., "Numerical recipes, the art of scientific computing," 
CDROM v2.10, Cambridge University Press, 2002. 


10 

American Institute of Aeronautics and Astronautics 



Table 1: Diagonal terms of linear stiffness obtained using five different modal bases. 
Shaded cells indicate bending eigenvalues. 



B 

BM 

CBC 

UBC 

UBC 9 (static) 

< 

1.318E+05 

1.318E+05 

1.318E+05 

1.318E+05 

1.318E+05 

< 

3.845E+06 

3.845E+06 

3.845E+06 

3.845E+06 

3.845E+06 

d] 

2.341E+07 

2.341E+07 

2.341E+07 

2.341E+07 

2.341E+07 

< 

8.079E+07 

8.079E+07 

8.080E+07 

8.079E+07 

8.079E+07 

< 

N/A 

4.988E+09 

N/A 

6.607E+00 

2.168E+02 

< 

N/A 

1.991E+10 

N/A 

1.029E+02 

8.083E+01 

d] 

N/A 

4.467E+10 

N/A 

3.208E+02 

7.326E+02 

< 

N/A 

7.908E+10 

N/A 

6.195E+02 

3.210E+02 


Table 2: Diagonal terms of cubic stiffness obtained using five different modal bases. 



B 

BM 

CBC 

UBC 

UBC 9 (static) 

i 1 2 

5.062E+02 

5.062E+02 

5.062E+02 

5.062E+02 

5.062E+02 

t 2 2 
^ 222^2 

3.615E+04 

3.615E+04 

3.615E+04 

3.615E+04 

3.615E+04 

T 3 2 

K<h 

2.488E+05 

2.488E+05 

2.488E+05 

2.488E+05 

2.488E+05 

T 4 2 

^444 ?4 

8.998E+05 

8.998E+05 

8.998E+05 

8.998E+05 

8.998E+05 

t 5 2 

N/A 

9.257E-05 

N/A 

-2.622E-08 

4.667E-07 

J 6 2 

^ 666^6 

N/A 

1.022E-03 

N/A 

5.324E-08 

5.568E-05 

T 7 2 

b m<h 

N/A 

-6.789E-04 

N/A 

-1.642E-07 

-5.252E-04 

T 8 2 

K*<h 

N/A 

4.139E-03 

N/A 

2.653E-07 

-4.691E-04 


Table 3: Selected off-diagonal terms of cubic stiffness obtained using five different modal bases. 



B 

BM 

CBC 

UBC 

UBC 9 (static) 

T 1 2 

Ka 

5.062E+02 

5.062E+02 

5.062E+02 

5.062E+02 

5.062E+02 


N/A 

1.270E-06 

N/A 

2.274E-1 1 

-8.560E-12 


N/A 

-1.173E+03 

N/A 

-2.801E-06 

-6.428E-07 

T 1 2 

b m 4, 

N/A 

2.592E-06 

N/A 

4.320E-14 

-1.638E-14 

r 5 2 

A 1U?1 

N/A 

-8.380E-02 

N/A 

1.730E-05 

2.610E-07 

*1 5 .5?1?5 

N/A 

-1.066E+03 

N/A 

-2.557E-06 

-7.970E-07 

b I m 

N/A 

4.680E-01 

N/A 

-5.137E-09 

9.984E-08 

i 5 2 

b ,A 

N/A 

9.257E-05 

N/A 

-2.622E-08 

4.667E-07 
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Figure 1: First anti-symmetric membrane and Figure 4: Fourth anti-symmetric membrane and 

companion modes. companion modes. 




igure 2: 


Second anti-symmetric membrane and 
companion modes. 


Figure 5: RMS transverse displacement as a function 
of distributed load. 




Figure 3: Third anti-symmetric membrane and 
companion modes. 


Figure 6: RMS membrane displacement as a 
function of distributed load. 
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Bending Stress PSD (psi 2 /Hz) Membrane Displacement PSD (in 2 /Hz) Transverse Displacement PSD (in 2 /Hz) 




Figure 7: Transverse displacement PSD for 
7.2xl0" 3 lb/in distributed load. 


Figure 10: Membrane stress PSD for 
7.2xl0" 3 lb/in distributed load 




Figure 8: Membrane displacement PSD for 
7.2xl0" 3 lb/in distributed load. 


Figure 11: Transverse displacement PSD for 
0.9216 lb/in distributed load. 
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Figure 9: Bending stress PSD for 
7.2xl0" 3 lb/in distributed load. 


Figure 12: Membrane displacement PSD for 
0.9216 lb/in distributed load. 
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RMS Transverse Displacement (in) Membrane Stress PSD (psi 2 /Hz) Bending Stress PSD (psi 2 /Hz) 




Figure 13: Bending stress PSD for 
0.9216 lb/in distributed load. 


Figure 16: RMS membrane displacement for 
thermally pre-buckled beam. 




Figure 14: Membrane stress PSD for 
0.9216 lb/in distributed load. 


Figure 17: Transverse disp. PSD for thermally pre- 
buckled beam (7.2xl0“ 3 lb/in distributed load). 




Figure 15: RMS transverse displacement for 
thermally pre-buckled beam. 


Figure 18: Membrane disp. PSD for thermally pre- 
buckled beam (7.2xl0 3 lb/in distributed load). 
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Figure 19: Bending stress PSD for thermally pre- 
buckled beam (7.2xl0' 3 lb/in distributed load). 
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Figure 22: Membrane disp. PSD for thermally pre- 
buckled beam (0.9216 lb/in distributed load). 



Figure 20: Membrane stress PSD for thermally pre- 
buckled beam (7.2xl0 3 lb/in distributed load). 



Figure 21: Transverse disp. PSD for thermally pre- 
buckled beam (0.9216 lb/in distributed load). 


Figure 23: Bending stress PSD for thermally pre- 
buckled beam (0.9216 lb/in distributed load). 



Figure 24: Membrane stress PSD for thermally pre- 
buckled beam (0.9216 lb/in distributed load). 
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Figure 25: Transverse displacement PSD for plate 
subjected to 106 dB acoustic loading. 


Figure 28: Membrane stress PSD for plate subjected 
to 106 dB acoustic loading. 




Frequency (Hz) 


Figure 26: Membrane displacement PSD for plate 
subjected to 106 dB acoustic loading. 


Figure 29: Transverse displacement PSD for plate 
subjected to 160 dB acoustic loading. 




Figure 27: Bending stress PSD for plate subjected to Figure 30: Membrane displacement PSD for plate 
106 dB acoustic loading. subjected to 160 dB acoustic loading. 
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Figure 31: Bending stress PSD for plate subjected to 
160 dB acoustic loading. 



Figure 32: Membrane stress PSD for plate subjected 
to 160 dB acoustic loading. 
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